#### weights are survey weights
plot_slopes_return_lmer_list <- function(
                        propensity="MARKET",
                        outcome="FAVOR>=3",
                        the_subset=sdat.under65.2,
                        the_subset2=NULL,
                        plot_title="Exchange Score\nAll respondents under 65",
                        pdf_name="pscore-by-survey-20190104-with-weights-lm.pdf",
                        y_axis="ACA Favorability",
                        .the_prefix=the_prefix,
                        color="purple",
                        color2=NULL, no_pdf=TRUE,
                        control_acafavor=FALSE,
                        no_polynomials = FALSE
                        ) {

    ## note - propensity is not the propensity score, it is training variable
    the_subset$propensity <- with(the_subset, eval(parse(text=propensity)))
    ## note - outcome is final outcome, not propensity outcome
    the_subset$outcome <- with(the_subset, eval(parse(text=outcome)))

    ## this can be used to apply insurance source scores to a holdout set
    ## but it ended up not being used in the paper
    if (!is.null(the_subset2)) {
        subset2_was_null <- FALSE
        the_subset2$propensity <- with(the_subset2, eval(parse(text=propensity)))
        the_subset2$outcome <- with(the_subset2, eval(parse(text=outcome)))
    } else {
        subset2_was_null <- TRUE
        the_subset2 <- the_subset
    }

    if (no_polynomials) {
    #### this is a model used to apply coefficients to small panel
    propensity_out_lm <- lm(
        propensity ~ EDUC + INCOME +
            AGEN +BLACK+HISP+ASIAN+MALE+RETIRED,
        data=subset(
            the_subset,
            DATN > 60 & !is.na(AGEN) & !is.na(INCOME) & !is.na(EDUC)
        ),
        weights=WEIGHT## WEIGHTS are survey weights
    )
    } else {
    #### this is the model used to produce the main results
    propensity_out_lm <- lm(
        propensity ~ poly(EDUC, 2) + poly(INCOME, 3) +
            poly(AGEN, 5)+BLACK+HISP+ASIAN+MALE+RETIRED,
        data=subset(
            the_subset,
            DATN > 60 & !is.na(AGEN) & !is.na(INCOME) & !is.na(EDUC)
        ),
        weights=WEIGHT## weights are survey weights
    )
    }
    ##
    #### this is a more interpretable model for the appendix
    propensity_out_lm_simple <- glm(
        propensity ~ scale(EDUC) + scale(INCOME) +
            scale(AGEN)+BLACK+HISP+ASIAN+MALE+RETIRED,
        data=subset(
            the_subset,
            DATN > 60 & !is.na(AGEN) & !is.na(INCOME) & !is.na(EDUC)
        ),
        family=poisson,
        weights=WEIGHT
    )

    ## produce propensity score
    the_subset2$PSCORE3 <- predict(
        propensity_out_lm,
        newdata=the_subset2,
        type="response"
    )
    propensity_sd <- sd(the_subset2$PSCORE3, na.rm=T)

    if (!control_acafavor) {
     lmer_out_all <- lmer(
        outcome ~ PSCORE3*I(PERIOD > 8.5) + (1|PERIOD),
        data=the_subset2,
        weights=WEIGHT
    )
     ##
    } else {
        lmer_out_all <- lmer(
        outcome ~ PSCORE3*I(PERIOD > 8.5) + I(FAVOR>=3) + (1|PERIOD),
        data=the_subset2,
        weights=WEIGHT
    )
    }

    unique_periods <- sort(
        unique(
            subset(
                the_subset2,
                !is.na(FAVOR) & !is.na(GOP) & !is.na(PSCORE3)
                )$PERIOD
        )
    )

    the_subset2$PSCORE3_orig <- the_subset2$PSCORE3

    the_subset2$PSCORE3 <- as.vector(
        scale(
            the_subset2$PSCORE3, scale=0.1 #for same scale across propensity scores
        )
        )

    if (!control_acafavor) {
     lmer_out_all_scaled <- lmer(
        outcome ~ PSCORE3*I(PERIOD > 8.5) + (1|PERIOD),
        data=the_subset2,
        weights=WEIGHT
    )
 } else {
        lmer_out_all_scaled <- lmer(
        outcome ~ PSCORE3*I(PERIOD > 8.5) + I(FAVOR>=3) + (1|PERIOD),
        data=the_subset2,
        weights=WEIGHT
        )
        }

    rmat <- matrix(NA,length(unique_periods),9)
    for(i in 1:(length(unique_periods))){
        dta.sub10 <- subset(
            the_subset2,
            PERIOD %in% unique_periods[i]
        )
        if (!control_acafavor) {
            lout <- lm(outcome ~ (PSCORE3), data=dta.sub10, weights=WEIGHT)## WEIGTHS are survey weights!
        } else {
            lout <- lm(outcome ~ (PSCORE3) + I(FAVOR>=3), data=dta.sub10, weights=WEIGHT)## WEIGTHS are survey weights!
        }
        rmat[i,1] <- unique_periods[i]
        rmat[i,2] <- dim(dta.sub10)[1]
        rmat[i,3] <- summary(lout)$coef[2,1]
        rmat[i,4] <- summary(lout)$coef[2,1] - 1.96*summary(lout)$coef[2,2]
        rmat[i,5] <- summary(lout)$coef[2,1] + 1.96*summary(lout)$coef[2,2]
        ##
        rmat[i,6] <- predict(
            lout, newdata=dta.sub10[1,] %>% mutate(PSCORE3 = 0), type="response"
        )
        rmat[i,7:8] <- predict(
            lout, newdata=dta.sub10[1,] %>% mutate(PSCORE3 = 0), type="response",
            interval="confidence"
        )[,2:3]
    }

    rmat[,9] <- as.numeric(rep(
        colnames(with(the_subset2, table(PERIOD, substr((DATE), 4, 5)))),
        apply(with(the_subset2, table(PERIOD, substr((DATE), 4, 5))), 2, function(x) sum(x>0))
    )[rownames(with(the_subset2, table(PERIOD, substr((DATE), 4, 5)))) %in% unique_periods])

    if (!no_pdf) {
    pdf(
        paste0(
            .the_prefix, "/figs/",
            pdf_name
        ),
        width=5, height=3.25
    )
    }
    par(mar=c(3,5,4,3))
    plot(x=rmat[,1]+0.5,y=rmat[,3],pch=16,
         xaxt="n",
         ylim=c(
             min(rmat[,4]),
             max(rmat[,5])
         ),
         xlab="",ylab="",type="n",
         bty="n",
         yaxt="n",
         xlim=c(0, 17)
         )
    coord.x <- c(9,16.75,16.75,9)
    coord.y <- c(-12,-12,12,12)
    polygon(coord.x,coord.y,col="darkgray", border=NA)
    par(new=T)
    plot(x=rmat[,1]+0.5,y=rmat[,3],pch=16,type="b",
         xaxt="n", ylim=c(
                       min(rmat[,4]),
                       max(rmat[,5])
                   ),
         xlab="Date",ylab=y_axis,cex.lab=1.3,
         main=plot_title,cex.main=1.25, bty="n", cex=1.5,
         ## yaxt="n",
         lwd=2,
         xlim=c(0, 17), cex.axis=1.2
         )
    ##
    ##
    axis(side=1, at=rmat[,1], labels=rmat[,9], cex=1.7, cex.axis=1.1)
    ##
    for(i in 1:dim(rmat)[1]){
        lb <- rmat[i,4]
        ub <- rmat[i,5]
        lines(x=c(rmat[i,1],rmat[i,1])+0.5,y=c(lb,ub), lwd=1.5)
    }
    mtext(side=3,at=5,"Pre-Exchanges", padj=0.25)
    mtext(side=3,at=13,"Post-Exchanges", padj=0.25)
    ##
    if (!is.null(color2)) {
        segments(
            y0=mean(rmat[rmat[,1]<9,3]),
            y1=mean(rmat[rmat[,1]<9,3]),
            x0=0, x1=9,
            col=color2, lwd=7, lty=1
        )
        segments(
            y0=mean(rmat[rmat[,1]>=9,3]),
            y1=mean(rmat[rmat[,1]>=9,3]),
            x0=9, x1=16.75,
            col=color2, lwd=7, lty=1
        )
    }
    segments(
        y0=mean(rmat[rmat[,1]<9,3]),
        y1=mean(rmat[rmat[,1]<9,3]),
        x0=0, x1=9,
        col=color, lwd=3, lty=1
    )
    segments(
        y0=mean(rmat[rmat[,1]<9,3]),
        y1=mean(rmat[rmat[,1]<9,3]),
        x0=9, x1=16.75,
        col=color, lwd=2, lty=3
    )
    segments(
        y0=mean(rmat[rmat[,1]>=9,3]),
        y1=mean(rmat[rmat[,1]>=9,3]),
        x0=9, x1=16.75,
        col=color, lwd=3, lty=1
    )
    if (mean(rmat[rmat[,1]>=9,3]) > mean(rmat[rmat[,1]<9,3])) {
        brackets(
            x1=16.75,
            x2=16.75,
            y1=mean(rmat[rmat[,1]>=9,3]),
            y2=mean(rmat[rmat[,1]<9,3]),
            col=color, xpd=F, lwd=3
        )
    } else {
        brackets(
            x1=16.75,
            x2=16.75,
            y1=mean(rmat[rmat[,1]<9,3]),
            y2=mean(rmat[rmat[,1]>=9,3]),
            col=color, xpd=F, lwd=3
        )
    }
    corners = par("usr")
    par(xpd = TRUE)
    text(
        x = corners[2]+1.1,
        y = mean(rmat[,3]),
        paste0(
            sprintf(
                "%.2f",
                round(mean(rmat[rmat[,1]>=9,3]) - mean(rmat[rmat[,1]<9,3]), 2)
            ),
            ifelse(
                abs(summary(lmer_out_all)$coefficients[4,3]) > 1.96,
            ifelse(
                abs(summary(lmer_out_all)$coefficients[4,3]) > 2.326,
                ifelse(
                    abs(summary(lmer_out_all)$coefficients[4,3]) > 3.09,
                    "***",
                    "**"
                    ),
                "*"
                ),
                ""
            )
           ),
        srt = 0, col=color,
        cex=1.25
    )
    ##
    if (!no_pdf) {
        dev.off()
        }

    return(
        list(
            lmer_out_all=lmer_out_all,
            lmer_out_all_scaled=lmer_out_all_scaled,
            propensity_out_lm=propensity_out_lm,
            propensity_out_lm_simple=propensity_out_lm_simple,
            propensity_sd=propensity_sd,
            rmat=rmat,
            out_data=the_subset2
        )
    )

}


plot_comparisons <- function(
                        plot_title="Uninsured vs Self-Purchased Score\nrespondents under 65, income over 40k",
                        pdf_name="selfinsure-vs-uninsure-pscore-by-survey-20190104-with-weights-lm.pdf",
                        y_axis="ACA Favorability",
                        .the_prefix=the_prefix,
                        color="purple4",
                        color1="orange4",
                        color2=gray(0.5),
                        mod_obj1,
                        mod_obj2
                        ) {


    pdf(
        paste0(
            .the_prefix, "/figs/",
            pdf_name
        ),
        width=5, height=3.25
    )

    rmat1 <- mod_obj1[["rmat"]]
    rmat2 <- mod_obj2[["rmat"]]

    par(mar=c(3,5,4,3))
    plot(x=rmat1[,1]+0.5,y=rmat1[,3],pch=16,
         xaxt="n",
         ylim=c(
             min(c(rmat1[,4], rmat2[,4])) - 0.02,
             max(c(rmat1[,5], rmat2[,5]))
         ),
         xlab="",ylab="",type="n",
         bty="n",
         yaxt="n",
         xlim=c(0, 17)
         )
    coord.x <- c(9,16.75,16.75,9)
    coord.y <- c(-12,-12,12,12)
    polygon(coord.x,coord.y,col=gray(0.8), border=NA)
    par(new=T)
    plot(x=rmat1[,1]+0.5,y=rmat1[,3],pch=21,type="b",
         xaxt="n", ylim=c(
             min(c(rmat1[,4], rmat2[,4])) - 0.02,
             max(c(rmat1[,5], rmat2[,5]))
                   ),
         xlab="Date",ylab=y_axis,cex.lab=1.3,
         main=plot_title,cex.main=1.25, bty="n", cex=1.5,
         lwd=2,
         xlim=c(0, 17), col=color1, bg=color
         )
    ##
    ##
    axis(side=1, at=rmat1[,1], labels=rmat1[,9], cex=1.5)
    ##
    for(i in 1:dim(rmat1)[1]){
        lb <- rmat1[i,4]
        ub <- rmat1[i,5]
        lines(x=c(rmat1[i,1],rmat1[i,1])+0.5,y=c(lb,ub), lwd=1.5, col=color)
    }
    mtext(side=3,at=5,"Pre-Exchanges", padj=0.25)
    mtext(side=3,at=13,"Post-Exchanges", padj=0.25)
    ##
    segments(
        y0=mean(rmat1[rmat1[,1]<9,3]),
        y1=mean(rmat1[rmat1[,1]<9,3]),
        x0=0, x1=9,
        col=color1, lwd=5, lty=1
    )
    segments(
        y0=mean(rmat1[rmat1[,1]<9,3]),
        y1=mean(rmat1[rmat1[,1]<9,3]),
        x0=9, x1=16.75,
        col=color1, lwd=4, lty=3
    )
    segments(
        y0=mean(rmat1[rmat1[,1]>=9,3]),
        y1=mean(rmat1[rmat1[,1]>=9,3]),
        x0=9, x1=16.75,
        col=color1, lwd=5, lty=1
    )
    ##
    if (mean(rmat1[rmat1[,1]>=9,3]) > mean(rmat1[rmat1[,1]<9,3])) {
        brackets(
            x1=16.75,
            x2=16.75,
            y1=mean(rmat1[rmat1[,1]>=9,3]),
            y2=mean(rmat1[rmat1[,1]<9,3]),
            col=color1, xpd=F, lwd=3
        )
    } else {
        brackets(
            x1=16.75,
            x2=16.75,
            y1=mean(rmat1[rmat1[,1]<9,3]),
            y2=mean(rmat1[rmat1[,1]>=9,3]),
            col=color1, xpd=F, lwd=3
        )
    }
    corners = par("usr")
    par(xpd = TRUE)
    text(
        x = corners[2]-5,
        y = mean(min(c(rmat1[1:8,3], rmat1[9:16,3])))-0.07,
        labels = "Self-insured",
        col = color1
        )
    text(
        x = corners[2]+1.1,
        y = mean(rmat1[9:16,3]),
        paste0(
            sprintf(
                "%.2f",
                round(mean(rmat1[rmat1[,1]>=9,3]) - mean(rmat1[rmat1[,1]<9,3]), 2)
            ),
            ifelse(
                abs(summary(mod_obj1[["lmer_out_all"]])$coefficients[4,3]) > 1.96,
            ifelse(
                abs(summary(mod_obj1[["lmer_out_all"]])$coefficients[4,3]) > 2.326,
                ifelse(
                    abs(summary(mod_obj1[["lmer_out_all"]])$coefficients[4,3]) > 3.09,
                    "***",
                    "**"
                    ),
                "*"
                ),
                ""
            )
           ),
        srt = 0, col=color1,
        cex=1.25
    )




    coord.x <- c(9,16.75,16.75,9)
    coord.y <- c(-12,-12,12,12)
    par(new=T)
    plot(x=rmat2[,1]+0.5,y=rmat2[,3],pch=21,type="b",
         xaxt="n", ylim=c(
                       min(rmat1[,4], rmat2[,4]) - 0.02,
                       max(rmat1[,5], rmat2[,5])
                   ),
         xlab="Date",ylab=y_axis,cex.lab=1.3,
         main=plot_title,cex.main=1.25, bty="n", cex=1.5,
         yaxt="n",
         lwd=2,
         xlim=c(0, 17), col=color2, bg=color
         )
    ##
    ##
    axis(side=1, at=rmat2[,1], labels=rmat2[,9], cex=1.5)
    ##
    for(i in 1:dim(rmat2)[1]){
        lb <- rmat2[i,4]
        ub <- rmat2[i,5]
        lines(x=c(rmat2[i,1],rmat2[i,1])+0.5,y=c(lb,ub), lwd=1.5, col=color)
    }
    mtext(side=3,at=5,"Pre-Exchanges", padj=0.25)
    mtext(side=3,at=13,"Post-Exchanges", padj=0.25)
    ##
    segments(
        y0=mean(rmat2[rmat2[,1]<9,3]),
        y1=mean(rmat2[rmat2[,1]<9,3]),
        x0=0, x1=9,
        col=color2, lwd=5, lty=1
    )
    segments(
        y0=mean(rmat2[rmat2[,1]<9,3]),
        y1=mean(rmat2[rmat2[,1]<9,3]),
        x0=9, x1=16.75,
        col=color2, lwd=4, lty=3
    )
    segments(
        y0=mean(rmat2[rmat2[,1]>=9,3]),
        y1=mean(rmat2[rmat2[,1]>=9,3]),
        x0=9, x1=16.75,
        col=color2, lwd=5, lty=1
    )
    if (mean(rmat2[rmat2[,1]>=9,3]) > mean(rmat2[rmat2[,1]<9,3])) {
        brackets(
            x1=16.75,
            x2=16.75,
            y1=mean(rmat2[rmat2[,1]>=9,3]),
            y2=mean(rmat2[rmat2[,1]<9,3]),
            col=color2, xpd=F, lwd=3
        )
    } else {
        brackets(
            x1=16.75,
            x2=16.75,
            y1=mean(rmat2[rmat2[,1]<9,3]),
            y2=mean(rmat2[rmat2[,1]>=9,3]),
            col=color2, xpd=F, lwd=3
        )
    }
    corners = par("usr")
    par(xpd = TRUE)
    text(
        x = corners[2]-5,
        y = mean(max(rmat2[1:8,3], rmat2[9:16,3]))+0.03,
        labels = "Uninsured",
        col = "black"
        )
    text(
        x = corners[2]+1.1,
        y = mean(rmat2[9:16,3]),
        paste0(
            sprintf(
                "%.2f",
                round(mean(rmat2[rmat2[,1]>=9,3]) - mean(rmat2[rmat2[,1]<9,3]), 2)
            ),
            ifelse(
                abs(summary(mod_obj2[["lmer_out_all"]])$coefficients[4,3]) > 1.96,
            ifelse(
                abs(summary(mod_obj2[["lmer_out_all"]])$coefficients[4,3]) > 2.326,
                ifelse(
                    abs(summary(mod_obj2[["lmer_out_all"]])$coefficients[4,3]) > 3.09,
                    "***",
                    "**"
                    ),
                "*"
                ),
                ""
            )
           ),
        srt = 0, col=color2,
        cex=1.25
    )
    ##


    dev.off()
}
